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Abstract 

Recent studies on the solvation of atomistic and nanoscale solutes indicate that a strong cou- 
pling exists between the hydrophobic, dispersion, and electrostatic contributions to the solvation 
free energy, a facet not considered in current implicit solvent models. We suggest a theoretical 
formalism which accounts for coupling by minimizing the Gibbs free energy of the solvent with 
respect to a solvent volume exclusion function. The resulting differential equation is similar to the 
Laplace- Young equation for the geometrical description of capillary interfaces, but is extended to 
microscopic scales by explicitly considering curvature corrections as well as dispersion and electro- 
static contributions. Unlike existing implicit solvent approaches, the solvent accessible surface is an 
output of our model. The presented formalism is illustrated on spherically or cylindrically symmet- 
rical systems of neutral or charged solutes on different length scales. The results are in agreement 
with computer simulations and, most importantly, demonstrate that our method captures the 
strong sensitivity of solvent expulsion and dewetting to the particular form of the solvent-solute 
interactions. 



I. INTRODUCTION 



Implicit solvent models are widely used in theoretical chemistry to study the solvation of 
biomolecular systems, as well described in the review of Roux and Simonsoni. They provide a 
more efficient, although generally less accurate, alternative to atomistically-resolved explicit 
solvent simulations. The solvation free energy in these models is usually split into nonpolar 
(np) and polar (p) terms, 

AG = AG np + AGp, (1) 

which are treated in separate energetic evaluations. The nonpolar term includes the energetic 
cost of cavity formation, solvent rearrangement, and solute-solvent dispersion interactions 
introduced when the uncharged solute is brought from vacuum into the solvent environment. 
The polar term describes the free energy of charging the mono- or multipolar solute in the 
dielectric medium. 

The nonpolar term is commonly approximated by surface area models, i.e. AG np ~ ^S, 
where 5* is the solvent accessible surface area£ and 7 is an energy per surface area constant, 
which is a priori not known but fit to atomistic simulations. The deficiencies of this simple 
surface area approach have been recognized and a further decomposition of the nonpolar 
term into cavity (cav) and van der Waals dispersion (vdW) terms has been proposed,-^ 
AG np = AG cav + AG v dw- This approach has shown improved results for the solvation of 
alkanes,^ the alanine peptide ® and nonpolar native and misfolded proteins. 7 The electro- 
static (polar) contribution of the solvation free energy is often approximated by generalized 
Born 8 (GB) or Poisson-Boltzmann- (PB) models. Both methods use a position-dependent 
dielectric constant,— assigned on the basis of the solute surface, which can be defined in 
several ways,-^ or defined implicitly by integration methods. It has been emphasized that 
all three contributions, AG cav , AG vdW , and AG p depend critically on the location of the 
solvent-solute interface. It has also been shown that the effective location of the solvent- 
solute interface can vary according to the local electrostatic^ and dispersion^ potentials. 
This suggests that interfacial, dispersion and electrostatic contributions should be coupled in 
implicit solvent approaches. The importance of capturing the right balance between nonpo- 
lar and electrostatic contributions in implicit solvation models was emphasized by Ashbaugh 
et al. in their study of amphiphiles.— 



The significance of nonpolar and polar coupling becomes even more evident when solva- 
tion is studied on length scales which are large compared to the solvent molecule (typically 
> 1 nm for water), where solvent dewetting ('drying') can occur. In this mechanism, first 
envisioned by Stillinger,— the solvent molecules tend to move away from the surface of 
a large nonpolar solute forming a liquid-gas like interface parallel to the solute interface. 
When the surfaces of two large solutes come together dewetting can be amplified due to the 
gain of interfacial free energy (by decreasing the total liquid-vapor interface area) giving rise 
to a strong effective attraction . 15 ' 1 ^ 1 ? Early evidence of confinement-induced dewetting was 
given only by explicit water simulations for smooth plate-like solutes with a purely repulsive 
solute-solvent interaction.^ More recently, however, it has been demonstrated in varying 
degrees in several systems with attractive solute-solvent interactions including smooth par- 
allel plate-like solutes,— atomistically resolved paraffin-plates,— graphite-plates 21 , carbon 
nanotubes,— and hydrophobic ion channels . 23124125 

Several of these studies indicated that the magnitude of dewetting is sensitive to the na- 
ture of the solute-solvent attractive dispersion interactions . 1 ? i2( - >i21 A similar sensitivity was 
found in systems where the solutes carry charges or are exposed to an external electric field, 
e.g. electrostatic interactions have been shown to strongly affect the dewetting behavior of 
hydrophobic channels 2 ^ 2 ^ 2 ^ and hydrophobic spherical nanosolutes . 29l3 ° Furthermore, two 
recent simulations of proteins supported the importance of solvent dewetting and its sensi- 
tivity in realistic biomolecular systems. First, a simulation of the two-domain BphC enzyme 
showed that the region between the two domains was completely dewetted when vdW and 
electrostatic interactions were turned off, but accommodated 30% of the bulk density with 
the addition of vdW attraction (water was found mainly at the edges of the considered vol- 
ume, while the central region was still empty), and 85-90% with the addition of electrostatic 
interactions.^ Second, Liu et al. observed a clear dewetting transition in the simulation of 
the collapse of the melittin tetramer, which was strongly sensitive to the type and location 
of the hydrophobic residues around the dewetted region? 3 ^ 

Considering the aforementioned studies, we postulate that coupling of the nonpolar and 
polar solvation contributions in implicit solvent models is crucial for an accurate determi- 
nation of solvation free energies without too many system-dependent fit parameters. We 
suggest a general theoretical formalism in which the particular energetic contributions are 
coupled. Similar to the approach of Parker et al. in their study of bubble formation at hy- 



drophobic surfaces,— we express the Gibbs free energy as a functional of the solvent volume 
exclusion function,^ and obtain the optimal solute surface via minimization. As we will 
show, this minimization leads to an expression which is similar to the Laplace- Young equa- 
tion for the description of macroscopic capillary interfaces,- 3 ^ but is generalized to explicitly 
include curvature corrections and solvent-solute interactions, i.e. short-range repulsion (ex- 
cluded volume), dispersion, and electrostatics. This extension of the Laplace- Young theory 
allows a geometric description of solvation on mesoscopic and microscopic scales. Related 
approaches in other fields are the Helfrich description of vesicle and membrane surfaces , ?9i?7 
wetting in colloids and granular media , 3 ^ 38 and functional treatments of electrowetting.^ 

While most implicit solvent approaches define the solute surface with a geometrical evalu- 
ation of the molecular surface, vdW surface, or canonical solvent accessible surface (SAS)r*^ 
it is an output of our theory. The surface obtained by minimizing our free energy functional 
will, in general, be very different than the aforementioned established surface definitions. In 
particular, our solvent accessible surface should not be confused with the canonical SAS, 2 
which is simply the envelope surrounding probe-inflated spheres. Similarly, phenomeno- 
logical continuum theories applied to solvent dewetting always assume a certain, simpli- 
fied geometry for the dry region, e.g. a cylindrical volume for system like hydrophobic ion 
channels , 24 i 28 i 40 plate-like particles,—^ 9 or two hydrophobic spherical solutes.— For a few sim- 
ple systems this might be a valid approximation but for more complicated solute geometries 
the shape of the dewetted volume is unknown and a different approach, as suggested in this 
work, is necessary. We expect our formalism to be particularly useful in solvation studies of 
large protein assemblies where the hydrophobic surfaces are highly irregular and laced with 
hydrophilic units , 42 > 43 and for which a unified description of hydration on different length 
scales is important.— Another potential application is the solvation of super hydrophobic 
nanosolutes^ 4 - and wetting/dewetting in near-critical colloidal mixtures. 38 

A brief summary of our work has been published elsewhere.^ Here we present more 
challenging test cases and an expanded discussion of the approximations and limitations 
of this model. The rest of the paper is organized as follows: In section II we present our 
theoretical formalism and chosen approximations. In section III we first verify that our 
method can describe solvation on molecular scales with noble gases, ions, and small alkanes. 
We then demonstrate that it captures the strong sensitivity of dewetting and hydrophobic 
hydration to specific solute-solvent interactions on larger scales with two alkane-assembled 



spheres. In section IV we conclude with some final remarks. 



II. THEORY 

A. Basic formalism 

Let us consider an assembly of solutes with arbitrary shape and composition surrounded 
by a dielectric solvent in a volume W. Furthermore, we define a subvolume (or cavity) V 
empty of solvent for which we can assign a volume exclusion function given by 



We assume that the surface surrounding the volume is continuous and closed, i.e. has no 
boundary. The absolute volume V and surface area S of V can then be expressed as func- 
tional of v(r) via 



where V = Vf is the usual gradient operator with respect to the position vector f and 
|Vf(r)| gives a 5-function-like contribution only at the volume boundary. The expression 
d 3 r |Vf(r)| = dS can thus be identified as the infinitesimal surface element. In this contin- 
uum solvent model the solvent density distribution is simply p{r) = p v(f), where po is the 
bulk density of the solvent at the desired temperature and pressure. Local inhomogeneities 
of the solvent density, apart from the zero to p transition at the volume boundaries, are 
neglected. The solutes' positions and conformations are fixed such that the solutes can be 
considered as an external potential to the solvent without any degrees of freedom. 

As motivated before, we suggest expressing the Gibbs free energy G[v] as a functional of 
the volume exclusion function v(r), and obtaining the optimal solute volume via minimiza- 
tion 




1 else. 



for f G V; 



(2) 




(3) 



5G[v]/Sv(r) = 0, 



(4) 



where S../Sv denotes the functional derivative with respect to the function v. We adopt 
following ansatz for the Gibbs free energy of the solvent: 



Let us discuss each of the terms in Eq. (jnj) in turn. The first term, G pr [i>], proportional to 
the volume V, is the energy of creating a cavity in the solvent against the difference in bulk 
pressure between the liquid and vapor phase, P = Pi — P v . For water in ambient conditions, 
which is close to the liquid-vapor transition, this term is relatively small and can generally 
be neglected for solutes on molecular scales. The second term Gmt[f] describes the energetic 
cost due to solvent rearrangement around the cavity interface with area S in terms of a free 
energy /surface area functional j(f; [v]). This interfacial energy penalty is thought to be the 
main driving force behind hydrophobic phenomena.— 7 is a solvent specific quantity that 
also depends on the particular topology of the cavity-solvent interface, i.e. it varies locally 
in space and is a functional of the volume exclusion function 7 = 7(r; [v]).— The exact form 
of this functional is not known. 

The third term, G ne [w], is the total energy of the non- electrostatic solute-solvent inter- 
action given a solvent density distribution pov(r). The potential U(r) = ^2nUi{f — r*j) is 
the sum of the (short-ranged) repulsive exclusion and (long-ranged) attractive dispersion 
interaction between each solute atom i at position r*j and a solvent molecule at r. Classical 
solvation studies typically represent the interaction Ui as an isotropic Lennard- Jones (LJ) 
potential, 



with an energy scale e, length scale ex, and center-to-center distance r. Using the form of 
(JHJ) implies that v(f) is defined with respect to the LJ-centers of the solvent molecules. 

The last term, G e sH, describes the total energy of the electrostatic field and the mobile 
ions in the system expressed by the local electrostatic potential \I/(r) assuming linear response 
of the dielectric solvent, 47 Similar to j(f] [v]), the position-dependent dielectric constant 



G[v] = G pr [v] + G int [v] + G ne [v] + G es [v] 



(5) 





(6) 



e(f) = e(f; [v]) depends on the geometry of v(r) with an unknown functional form. A(r) is 
the fixed charge density distribution of the solutes and the local energy density of the mobile 
ions is*^ 

t/ mi (r) = fc B T^{exp[-/%*(f)] - 1}. (7) 

3 

with the thermal energy k^T = (3~ l . Variation of (0) for a fixed v(f) with respect to \I/(r) 
yields the Poisson-Boltzmann equation*^ 

PB(f) = = V • [e e(f; M)W(f)] + A(f) 

+ ^(rl^^-exph^lrl, (8) 

i 

where gj and pj are the charge and concentration of the mobile ion species j. Note that 
the ionic charge density in (JBJ) is multiplied by v(f) to account for the fact that ions usually 
cannot penetrate the volume empty of polar solvent due to a huge free energy penalty. We 
remark that the treatment of the electrostatics in our theory has the same limitations as 
other implicit models using PB, for instance when describing highly charged or strongly 
correlated electrolyte systems. In contrast to PB/SA models however, the dielectric bound- 
ary is optimized such that it responds to the local nonpolar and polar potential; it is not 
assumed beforehand. 

Let v m i n (f) be the exclusion function which minimizes the functional (jSJ). Then, the 
resulting Gibbs free energy of the system is given by G[t> m in]- The solvation free energy AG 
is the reversible work to solvate the solute and is given by 

AG = G[v min ) - G , (9) 

where Go is a constant reference energy which can refer to the pure solvent state and an 
unsolvated solute. The potential of mean force (pmf) along a given reaction coordinate 
x (e.g. the distance between two solutes' centers of mass) is given, within a constant, by 
G[f m in], where v m i n (r) must be evaluated for every x. In order to proceed we will need valid 
approximations for 7(r; [v]) and e(f; [v]) with which f m i n (^*) can be calculated by explicitly 
minimizing our free energy functional (0) according to 



B. Approximations for y(r; [v]) and e(r; [v]) 

Let us start with a possible description of 7(r*; [v]). For a planar macroscopic interface 
the parameter 7 is usually identified by the surface tension of the solvent adjacent to the 
second medium. This surface tension obviously depends on the microscopic interactions 
between the medium and the solvent and is generally decreased by attractive dispersion or 
electrostatic contributions. It seems that microscopic interactions are adequately represented 
by a macroscopic quantity like 7 if their range is much smaller than the investigated length 
scales, such as the radii of curvature or mean particle distances. The effect of the microscopic 
interactions are then absorbed in 7. This has been exemplified with free energy estimates 
for the solvation of large, neutral plate-like or spherical alkane-assembled solutes , 12 ' 1 ? For 
the description of solvation on smaller length scales, however, it seems important to separate 
the free energy into a part which accounts for the formation of a cavity and a part which 
describes the dispersion interactions explicitly. 3 Furthermore, it has been shown that the 
water liquid-vapor surface tension 7^ is the asymptotic value of the solvation free energy 
per surface area for hard spherical cavities in water in the limit of large radii . 1 ^ 4 ? These 
considerations motivate our choice of the second and third term in the functional (jHJ) and 
lead to the assumption 7 = 7i v in the limit of vanishing curvatures. 

The surfaces of realistic (bio)molecules, however, display highly curved shapes, so j(f; [v]) 
will strongly depend on the interface geometry around r in a complicated fashion. In the 
following we make a local curvature approximation, i.e. we assume that j(f; [v]) can be 
expressed solely as a function of the local mean curvature 

H(f) = (« x (r) + K 2 (r)) /2 = R{f)~\ (10) 

where R(r) is the radius of mean curvature and K\{r) and K%(f) are the local principal 
curvatures of the interface.— The mean curvature H is only defined at the boundary of v(r). 
We have chosen the convention in which the curvatures are positive for convex surfaces (e.g. 
a spherical cavity) and negative for concave surfaces (e.g. a spherical droplet). 

The curvature dependence of the liquid-vapor surface tension has been a long standing 
subject of research and is still under steady discussion . 3 ^ 1 ^ 2 For water, which is close to 
the critical point under ambient conditions, 7 is argued to be nonanalytical in curvature.— 
The first order correction term, however, is likely to be linear in curvature as predicted by 
scaled-particle theory— the commonly used ansatz to study the solvation of hard spherical 



cavities. Although this result is only strictly valid for the case of spherical particles, we 
assume that it can be applied to local mean curvatures such that 7(f; [v]) reduces to the 
function 



where 5 is the Tolman length, which is expected to be of molecular size.— In our study 
we assume 5 is constant and positive, while the curvature can be positive or negative as 
defined above. Note that this leads to an increase of surface tension for concave surfaces 
in agreement with the geometrical arguments of Honig et a/.— in their solvation study of 
alkanes. It has been shown by computer simulations of growing a hard spherical cavity in 
water that (|TT|) predicts the interfacial energy rather well for radii > 3A.— A major drawback 
of approximation (JTTJ) is that it gives unphysical results if the radius of mean curvature is 
smaller than twice the Tolman length, \R\ < 25. It yields negative and diverging surface 
tensions for convex and concave surfaces, respectively. The latter is not possible due to the 
finite size of the solvent molecules. Thus, care has to be taken with approximation (JTTJl 
when investigating systems which can exhibit radii of curvature \R\ < 25. 

Let us now turn to electrostatics. The most common approximation for the position- 
dependent dielectric constant is proportional to the volume exclusion function f(f),- such 
that the functional e(f; [v]) reduces to, 



where e v and ei are the dielectric constants inside and outside the volume V, respectively. 
Eq. (|12jl is valid only in the limit of large solute sizes when the molecular size of the solvent 
is negligible. For charged solutes on a molecular scale, let's say mono- or polyvalent ions, 
two difficulties arise. First the electric field close to the highly curved solutes can be strong 
enough for the dielectric constant to be field dependent. This formally affects the form of 
the electrostatic term in the free energy functional which assumes a linear response of the 
solvent. An improvement for continuum models along these lines has been proposed by Luo 
et alM Second, the effective position of the dielectric boundary is known to depend on the 
sign of the solute charge for asymmetric solvent molecules like water. This expresses itself, 
for instance, in different Born radii for two equally charged ions which have exact the same 
LJ parameters but a different sign of charge. A reasonable improvement of (|T2*|) would be to 



7 (f) = 7iv(l - 25H(f)) 



(11) 



e(f) = e„ + v(r)(ei - e v ), 



(12) 



shift the dielectric boundary at r parallel to the volume boundary by a potential dependent 
amount Ar = '(f))n(r): 

e(f) =e„+v(r- Ar)(e,-e„), (13) 

where n is the unit normal vector to the interface. We do not attempt however, to find an 
approximation for the function in this work and postpone this investigation to later 
studies. For further illustration of our approach we content ourselves with the approxima- 
tions (fTTj) for 7(r; [v]) and (fT2j) for e(r; [v]). 



C. Minimization of the free energy functional 

For the functional derivative of the interfacial term, GmtM, we utilize 

— [ d 3 r \Vv(r)\ 

and 

■f- [ d 3 rH \ Vv(r)\ 
ov J w 

which has been derived in detail by Zhong-can and Helfrich by means of differential 
geometry— The variable K(r) = Ki(r)«2(f) is the Gaussian curvature of the interface, 
which is an intrinsic geometric property of v. Plugging in approximations (jllj) and (fT2j) into 
(jSJ, and minimizing with (0J), we obtain, 

= de(f) = P + 2 7 i v [H(r) - 6K(r)} - p U(r) 

- |[V*(f)c(f)] 2 f- - -) -U mi (r). 

(16) 

Eq. (fTB*j) is a partial second order differential equation (de) for the optimal exclusion func- 
tion v m i n (f) expressed in terms of pressure, curvatures, short-range repulsion, dispersion, and 
electrostatic terms, all of which have dimensions of energy density. It can also be interpreted 
as a mechanical balance between the forces per surface area generated by each of the par- 
ticular contributions. Thus, in our approach the surface shape and geometry, expressed by 
H and K, are directly related to the inhomogeneous potential contributions. The constant 
solute charge density X(r) does not appear explicitly in (fTB^) but is implicitly considered in 



8v 



dS = -2H 



(14) 



. , dS H = -K, 



(15) 



the PB equation (jHJ), which must be solved simultaneously. If curvature correction (if -term) 
and the last three energetic terms are neglected one obtains the Laplace- Young equation, 

P = -2 llv H, (17) 

which is exclusively used for the shape description of macroscopic capillary and interfacial 
phenomena in conjunction with appropriate boundary conditions, e.g. prescribed liquid-solid 
contact angles at the solid surfaces.— In our description the boundary conditions are provided 
by the constraints given by the short-ranged repulsive term in U (f) , and the distribution of 
dispersion and electrostatics, allowing an extrapolation of the Laplace- Young description to 
mesoscopic and microscopic scales. Notice that in our approach the solvent is treated as a 
continuum while the solute is explicitly resolved. One could use a coarse-grained treatment 
for the solute by including the appropriate non-electrostatic and electrostatic interactions 
in ©. 

The solution of (fTBj) requires an appropriate parametrization, i.e. coordinate representa- 
tion, for the curvatures H and K, such that the equation is expressed as a function of the 
vector r*and its first and second derivatives in space. Analytical solutions to the much sim- 
pler Eq. (fTTj) and thus to (fTfij) are only available for systems with very simple geometries.^ 
Thus we use numerical solutions of (j!6|) in the following to further illustrate our theory. 

III. APPLICATIONS 

First, we will consider the solvation of microscopic solutes such as noble gases, simple 
alkanes, and ions which can be treated as neutral or charged Lennard- Jones spheres. Then, 
we will investigate alkane assemblies on a larger scale where interfacial and dewetting effects 
are much more dominant. For simplicity and a better transparency of the results, mobile 
ions will be neglected in these illustrations. 

A. One Lennard-Jones sphere 

In this section we compare our approach to results from SPC and SPC/E explicit sol- 
vent simulations.-^ We refrain from comparing to real experiments since approximations in 
computer experiments are more easily controlled and the LJ parameters of the solutes are 
commonly parametrized to yield accurate results in classical simulations. 



For a spherical solute with a charge Q homogeneously distributed over its surface, the 
functional (jSJ) with approximations (fTTj) and (|T2|) and no mobile ions reduces to a function 
of R, the radius of the sphere empty of solvent. The solvation free energy is 



AG{R) = AG pr {R) + AGint(i2) + AG ne (R) + AG es {R) 
= -irR 3 P + 4irR 2 llv (l- — 

Anr 2 dr p U L j(r) 

Note that the last term in (JTSj) is equivalent to the Born electrostatic solvation free energy.— 
Recently, Manjari et al. have presented a very similar expression for the solvation of a 
charged spherical cavity on the basis of a minimization principle and have investigated the 
variation of R with thermodynamic conditions. 58 Differentiation of (|18|) with respect to R 
and subsequent division by 4ttR 2 yields 

which is in accord with Eq. given sphere-like curvatures, H = 1/R and K = 1/R 2 . We 
can now calculate the solvation free energies of simple spherical solutes, such as noble gases 
or ions. The free parameters in Eq. (j!9|) are the pressure P, Tolman length 5, liquid-vapor 
surface tension 7i v , and dielectric constants e v and t\. 



I. One neutral LJ sphere 

First, let us focus on uncharged spheres, for which the electrostatic term in (fT5|) can 
be neglected. We compare the results from our theory to those calculated by Hummer et 
al^ for neutral LJ spheres in SPC water, and those calculated by Paschek^ for noble gases 
in SPC and SPC/E water. The solute- water LJ parameters a and e are summarized in 
Tab. I. The surface tension 7i v was set to that estimated for SPC and SPC/E water at 
300K, 7 lv = 65mJ/m 2 and 7i v = 72mJ/m 2 , respectively.—^ The pressure is fixed to latm. 
Finally, the remaining free parameter 5 was fit to reproduce the simulation solvation free 



energies exactly. The solvation free energies from simulation AG s i m and best fit Tolman 
lengths 5m are shown in Tables I and II for the SPC and SPC/E models, respectively. 

Before we discuss the results, let us compare the particular energy contributions AGi(R) 
with i = pr, int, ne for Na° (plotted in Fig. QJ. As anticipated, the pressure term AG pr (R) 
with P = latm is negligible compared to the other contributions. The interfacial term 
AG int (i?) increases with the cavity radius R. The integrated LJ-interaction term AG ne (R) 
shows long-range attraction and a steep short-ranged repulsion with a minimum at R = 
a(Na°) = 2.85A. The total solvation free energy for the Na° shows a single minimum at 
R min = 2.32A with AG = 9.2kJ/mol for a 5 bf =0.79A. 

The results for the other LJ-spheres, summarized in Tab. I and II, reveal several note- 
worthy observations. First, the best fit Tolman lengths 5m range from 0.76 A to 1.00 A ; 
they are not only of molecular size, as expected, but are approximately half the LJ-radius 
of a SPC or SPC/E water molecule. Second, the 5m values for noble gases in SPC/E wa- 
ter (Tab. II) are approximately 10% larger than those in SPC water (Tab. I). This is in 
qualitative agreement with the results of Huang et al. who measured 5=0.76 ± 0.05A and 
5=0.90 ± 0.03A for SPC and SPC/E, respectively, by fitting Eq. dHJ) to the hydration free 
energy of hard spheres with varying radii.— 

Third, the quite accurate data of Paschek demonstrate a systematic increase of 5m with 
solute size. The inability of our theory to be fit by one fixed constant 5m points to the 
anticipated fact that Eq. (fTT]) can not capture strong curvature effects accurately and will 
have to be refined for small solutes. Despite this shortcoming, these results show surprisingly 
good agreement; if we assume a fixed delta, for instance 5 = 0.9lA for all noble gases in the 
SPC data of Paschek, our theory predicts results within 15% of the simulation data. Finally, 
we observe that the effective optimal sphere radius -R m i n is always smaller than the radius 
of the canonical SAS with a typical probe radius of -R m i n < (cr ss /2 + 1.4 A) ~ a, but 

larger than the vdW surface, -R m i n > cr ss /2, where a ss is the solute-solute LJ-length.^ 

2. One charged LJ sphere 

Let us now turn to charged Lennard- Jones spheres (ions) also examined in the paper by 
Hummer et al. with SPC water simulations. We assume 5 to be the fixed by the previously 
obtained <5bf values for uncharged spheres. The dielectric constants are set to e v — 1 and q = 



65, in accord with SPC water.— The electrostatic contribution AG es (R) and the total AG(R) 
for Na + are shown in the inset of Fig. The electrostatic contribution decreases the optimal 
radius to -R m i n = 1.83 A giving a solvation free energy of AG = — 334kJ/mol. In fact, the 
optimal sphere radius -R m i n is always considerably smaller for the charged solutes (Tab. Ill) 
than for their neutral counterparts (Tab. I). This is caused by the strong compressing force 
of the polar solvent attempting to penetrate the low dielectric cavity. The solvation free 
energies from theory AG and those from simulation AG s ; m are also shown in Tab. III. While 
our theory describes the hydration free energies for positively charged ions within 15%, 
it considerably underestimates those of the negative ions. This qualitative disagreement 
between positive and negative ions was expected since the Born radii for anions are always 
smaller than those for cations, a consequence of the different solvation structure around 
charged solutes with opposite signs. As mentioned in the previous section, the position of the 
dielectric boundary has to be refined for accurate estimates of the electrostatic contribution 
to the hydration free energy. If we apply the correction (jTBj) to the dielectric boundary 
with a simple, potential-independent shift £ + = — 0.25A for positive and = — 1.05A for 
negative spheres such that the dielectric boundary has a radius of R + £± < R, improved 
values (AGg in Tab. Ill) are obtained which reproduce all simulation values within 10%! 

B. Linear alkanes 

Let us now consider simple polyatomic molecules, such as ethane, propane, or butane in 
a cylindrically symmetric one- dimensional (ID) chain conformation. Other conformations 
will be neglected. The symmetry of these systems allows us to express the volume exclusion 
function v(r) of the enveloping surface by a one dimensional shape function r(z), where z is 
the coordinate on the symmetry axis and r the radial distance to it. The full surface in three- 
dimensional space is obtained by revolving the shape function r(z) around the symmetry 
axis. Technical details are given in the Appendix. 

The LJ parameters for ethane and methane are the same as those used by Ashbaugh et 
alJ& in their SPC simulation of linear alkanes (see Tab. I). The simulation solvation energy 
of the spherical methane, AG = 10.96kJ/mol, can be reproduced with 5m = 0.85A. Solving 
the cylindrically symmetric problem for ethane using the same 5, we obtain a fit-parameter 
free AG = 11.40kJ/mol, which is only 7% larger than the simulation results. Alternatively, 



5bi = 0.87A reproduces the simulation energy exactly. This is surprisingly good agreement 
considering the crudeness of our curvature correction and the fact that the large curvature 
of the system varies locally in space. The curvature and shape functions are plotted in 
Fig. El together with the vdW surface and the canonical S AS obtained from rolling a probe 
sphere with the typically chosen radius r p = 1.4A over the vdW surface.— Away from 
the center of mass \z\ > lA the curvatures follow the expected trends for the spherical 
surfaces: H — 1/R and K = 1/R 2 with R ~ 3.1 A . The optimal surface resulting from our 
theory is smaller than the canonical SAS and smooth at the center of mass (z = 0) where 
the canonical SAS has a kink. Thus our surface has a smaller mean curvature at z = 
and an almost zero Gaussian curvature, which is typical for a cylinder geometry with one 
principal curvature equal to zero. These results may justify the use of smooth surfaces in 
coarse-grained models of closely-packed hydrocarbons, a possibility we will explore in the 
following section with solvation on larger length scales where dewetting effects can occur. If 
we repeat the above calculation for propane and butane (three and four LJ-spheres, see also 
Tab. I for parameters) we need 5m = 0.94A and 5m = 0.96A, respectively, to reproduce the 
simulation results exactly. The increasing difference in 5m compared to methane and ethane 
is likely due to contributions from other than cylindrically symmetric conformations which 
were ignored in our analysis. 

C. Two spherical nanosolutes 

1. Model 

Let us now consider two spherical solutes which represent homogeneously assembled CH2 
groups with a uniform density p=0.024A~ 3 up to a radius Rq = 15A, defined by the maximal 
distance between a CH2 center and the center of the solute. Integration of the CH2-water 
LJ interaction over the entire sphere yields a 9-3 like potential Ui(r) for the interaction 
between the center of the solute (i = 1,2) and a water molecule.— The intrinsic, nonelec- 
trostatic solute-solute interaction U ss can be obtained in a similar fashion. The CH 2 -water 
LJ parameters, e = 0.5665kJ/mol and a = 3.52A, are taken from the OPLSUA force-field^ 
and are similar to those used by Huang et al. in their study on dewetting between paraffin 
plates.-^ Minimizing Eq. (|T%|) for just one sphere we obtain an optimal solvent excluded 



radius of -R m i n ~ 17.4A, which is Ro + 2.4A. Since we are also interested in the effects of 
electrostatic interactions we place opposite charges ±Ze, where e is the elementary charge, 
in the center or on the edge of the two spheres. Poisson's equation is simultaneously solved 
on a two-dimensional grid in cylindrical coordinates. Numerical details are given in the 
Appendix. 

The solvation of the two solutes is studied for a fixed surface-to-surface distance which 
we define as s = r 12 — 2R Q , where r 12 is the solute center-to-center distance. The effective 
surface-to-surface distance defined by the accessibility of the solvent centers is thus s ~ 
ri2 — 2R mm = sq — 4.8A. In the following we focus on a separation distance of so = 8A to 
investigate the influence of different energetic contributions on the shape function, r(z), 
and the curvatures, K(z) and H(z). For s = 8 A, it follows that s ~ 3.2A, such that 
two water molecules could fit between the solutes on the z-axis. We systematically change 
the solute-solute and solute-solvent interactions, as summarized in Tab. I. We begin with 
only the repulsive part of the nonelectrostatic interaction £/j(r) in system I, and then add a 
curvature correction with 5 = 0.75A, vdW attractions, and sphere-centered charges Z = 4 
and Z = 5 in systems II-V, respectively. To study the influence of charge location, we 
reduce the magnitude of each charge in system VT to Z — 1 and move them to the edge 
of the spheres on the symmetry axis such that they are 8A apart (indicated by arrows in 
Fig. El The surface tension and dielectric constants of the vapor (solute) and liquid are fixed 
to 7i v = 72mJ/m 2 , e v = 1, and q = 78, respectively. 

2. Behavior of the shape function 

The results for the curvatures and shape function, defined by r(z), for systems I- VI are 
shown in Fig. El Away from the center of mass (\z\ > 10 A), systems I- VI show very little 
difference. The curvatures are H = l/R and K = 1/R 2 with R ~ 17.4A. Close to the 
center of mass ~ 0), however, the influence of changing the parameters is considerable. 
In system I, Eq. (|T^|) reduces to the minimum surface equation H(z) = for z ~ 0. For two 
adjacent spheres the solution of this equation is the catenoid r(z) ~ cosh(z), which features 
zero mean curvature (ki and n 2 cancel each other) and negative Gaussian curvature. As a 
consequence, the system exhibits a vapor bubble bridging the solutes, i.e. water is removed 
from the region between the spheres even though it fits there. This dewetting is driven by 



the interfacial term G- m t which always favors minimizing the liquid-vapor interface. 

When curvature correction is applied (system II), the mean curvature becomes nonzero 
and negative (concave) at z ~ 0, while the Gaussian curvature grows slightly more negative. 
Thus the total enveloping surface area becomes larger and the solvent inaccessible volume 
shrinks, i.e. the value of the shape function at z ~ decreases. Turning on solute-solvent 
dispersion attraction amplifies this trend significantly as demonstrated by system III. Mean 
and Gaussian curvatures increase fivefold, showing strongly enhanced concavity, and the 
volume empty of water decreases considerably, expressed by r(z = 0) ~ 10. 7A dropping to 
r(z = 0) ~ 6.3A. These trends continue with the addition of electrostatics in system IV. 
When the sphere charges are further increased from Z = 4 to Z = 5 (system IV— >V), we 
observe a wetting transition: the bubble ruptures and the shape function jumps to the solu- 
tion for two isolated solutes, where r(z ~ 0) = 0. The same holds when going from III to VI, 
when only one charge unit, Z = 1, is placed at each of the solutes' surfaces. Importantly, this 
demonstrates that the present formalism captures the sensitivity of dewetting phenomena to 
specific solvent-solute interactions as reported in previous studies . 2 ' ' 112112 ^ 2 ? 12 ^ 1 ^ 2 Note that 
the optimal shape function at \z\ ~ ±2 A is closer to the solutes in VI compared to V due 
to the proximity of the charge to the interface. Clearly, the observed effects, particularly the 
transition from III to VI, cannot be described by existing solvation models which use the sur- 
face area (GB/SA or PB/SA) 1 or effective surface tensions and macroscopic solvent-solute 
contact angle o?? 1 ^ as input. 

3. Potential of mean force 

The significant change of the shape function with the solute-solvent interaction has a 
strong impact on the potential of mean force (pmf) (or effective interaction) between the 
solutes 

W(s ) = G(s ) - G(oo) + U BB (s ). (20) 

Recall that U SB is the instrinsic dispersion interaction between the two solutes. Values of 
W(sq = 8A) are given in Tab. IV. From system I to VI the total attraction between the 
solutes decreases almost two orders of magnitude. Interestingly, the curvature correction 
(I — ^11) lowers W by a large 23.5fcsT, even though R ^> 5. The reason is that the mean radii 



of curvature between the spheres can assume values ~ 5, implying that curvature correction 
is also important for large solutes. A striking effect occurs when vdW contributions are 
introduced (II — >III) : the inter solute attraction decreases by approximately 28kBT while 
the dispersion solute-solute potential U ss (so = 8A) changes by only — 0.44fceT. Similarly, 
adding charges of Z = 5 (III — > V) at the solutes' centers or Z = 1 (III — > VI) at the 
solutes' surfaces decreases the total attraction by and 5k^T, respectively. Note that 

the total attraction decreases even though electrostatic attraction has been added between 
the solutes. The same trend has been observed recently in explicit water simulations of a 
similar system of charged hydrophobic nanosolutes.^Si 

Now we turn our attention to varying the intersolute distance. The pmfs and solute- 
solute mean forces F = dW(s )/ds versus a range of s are shown for systems I, II, III, 
and VI in Fig. HJ System I, with purely repulsive solute-solvent interactions, displays a 
strong attraction (W ~ — lSO&^T) at so = 2A which decreases, almost linearly, to zero at a 
distance so = 13. 5A where the system shows a wetting transition. The corresponding force 
is discontinuous at this critical distance. The steep repulsion at short intersolute distances 
(sq ~ 1.5A ) stems from the repulsive term of the LJ interaction between the solutes. Note 
that the intrinsic solute-solute interaction U ss (so), also shown in Fig.|H is almost two orders 
of magnitude smaller than the hydrophobic attraction. Adding the curvature correction 
in system II decreases the range and strength of the pmf by approximately 20%, which 
is significant and unexpected since Rq ^> 5. Adding dispersion attractions in system III 
decreases the range and strength of the hydrophobic attraction considerably, but it is still 
much stronger than the inter solute dispersion attraction U ss alone. When surface charges 
(Z = 1) are added in system VI, the range of hydrophobic attraction further decreases but 
the total attraction increases at short intersolute distances. This is due to the increasing 
size of the bridging bubble (r(z = 0) increases) as the two solutes approach each other, 
which decreases the high dielectric screening of the solute-solute electrostatic attraction. 
This again underlines the importance of coupling electrostatics and dewetting effects, as the 
electrostatic attraction (or repulsion) may be magnified by more than an order of magnitude 
when dewetting occurs. For charges with opposite sign this could be interpreted as the 
stabilization of a salt bridge due to dehydration.— Systems IV and V, not shown in Fig. HI 
exhibit the same qualitative behavior as system VI. 



4- Comparison of mean forces to previous MD simulations 



We continue considering the mean force between two nanosized solutes and compare our 
theory now to the MD simulations of Dzubiella et alM*& Their solute model slightly differs 
from the one used in the previous section: the solute-solvent interaction potential is purely 
repulsive and is given by Ui(r) = ksTir — Rq)~ 12 , while the solute-solute interaction is hard- 
sphere like with a hard sphere radius Rq. The solutes are neutral or carry opposite charges 
Q homogeneously distributed over the sphere volume. The simulations were carried out 
with the SPC/E model of water. In our theory, we fix the Tolman length to 5 = 0.90A as 
measured by Huang et al^ and the dielectric constant to q = 71 for SPC/E water The 
mean forces are shown in Fig. |5] for neutral spheres of radii i?o = 10 and 12A and oppo- 
sitely charged solutes with radius Rq = 10A and charge Q = 2e and 5e versus the solutes' 
surface-to-surface distance sq. Simulation and theory are in good, almost quantitative agree- 
ment, and show that our theory captures the decreasing range of the strongly hydrophobic 
attraction with decreasing radius and increasing charge due to suppressed dewetting. We 
emphasize that our theory is basically fit-parameter free for this system of large solutes. 
Fig. El also shows the theoretical mean force for the neutral Rq = 12 A solute using a smaller 
Tolman length 5 = 0.75A. The decrease of the Tolman length increases the depth and range 
of the solvent-mediated solute-solute attractive mean force by approximately 5%, showing 
a nonvanishing but only slight influence. 



IV. CONCLUSION AND FINAL REMARKS 



In summary, we have presented a novel implicit solvent model which couples polar and 
nonpolar solvation contributions by employing a variational formalism in which the Gibbs 
free energy of the system is expressed as a functional of the solvent volume exclusion function. 
Minimization of the free energy leads to a Laplace- Young like equation for the solvent 
excluded cavity around the solutes, which is extended to describe solvation on mesoscopic 
and microscopic scales. We have shown that the theory gives a reasonable description 
of the solvation of microscopic solutes, such as ions and alkanes. Improved accuracy will 
require further refinement of the curvature dependence of the surface tension 7(f; [v]) and 
the definition of the position-dependent dielectric constant e(r; [v]). Given the physically 



reasonable values of the parameters S and £ we found by fitting, we hope that extensions 
based on physical rational, e.g. given by complementary microscopic approache D 1 ^ 4 ' 4 ^? 1 ^ 
and further empirical corrections, will lead to an accurate fit-parameter free implicit solvent 
description. 

We have further demonstrated that on larger scales, where solvent dewetting can 
play an important role in solvation, our formalism captures the delicate balance be- 
tween hydrophobic, dispersive and electrostatic forces which has been observed in previous 
systems . 2 ^ 2112 ^ 2 ? 12 ? 1 ? 1 ^ 2 The dewetting in our model is driven by the interfacial term which 
favors minimizing the solute-solvent interface. A comment must be made here regarding the 
sensitivity of dewetting to the particular solvent-solute interactions. As recently argued by 
Chandler,— extended fluid interfaces near phase coexistence are often referred to as 'soft' 
because they can be deformed with only little or no free-energy change.— Our approach 
seems to account for this sensitivity since small changes of the constraints in the differential 
equation (jlfij) for the shape function, given e.g. by the dispersion potential close to the 
solute surface, can lead to a major deformation or even rupture (wetting transition) of the 
inter-solute, dewetted region. As we have shown, this can significantly change the pmf for 
the solutes. Thus we anticipate that slight changes in the geometry of a system, e.g. a slight 
concave or convex bending of two plate-like solutes , 20 ! 21 can lead to very different results for 
the dewetting magnitude and the pmf. 

The current illustrations utilized spherical and cylindrical symmetries. More complex 
molecules, such as proteins, will require solving the full three dimensional problem. Numer- 
ical algorithms for the calculation of interface evolution for more complicated geometries 
are provided by efficient level-set methods or fast marching methods.™ We believe that in 
the full three-dimensional (3D) case, our method will be much more efficient than other 
microscopic approaches which partly resolve the water structure and are able to describe 
dewetting effects, e.g. the Lum- Chandler- Weeks theory^ (LCW) or information theor y^?? 
(IT), as only a two-dimensional surface is sought rather than a 3D density distribution on a 
fine grid. We remark that LCW and IT do not consider electrostatic interactions and may 
benefit from our complementary approach. 
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Appendix A: Curvatures in cylindrical coordinates 

In our general parametrization for the shape function r(z) we express the radial coordinate 
by r = r(l) and the axial coordinate by z = z(l), as functions of the parameter I. Depending 
on the geometry of the considered system, / has to be conveniently chosen, for instance to 
be the arc length, or r(l) = I, or z(l) = I. In our illustration the most convenient choice is 
z(l) = I. The principal curvatures are generally given by^ 

-z' , s z'r"-z"r' 

Kl(r ' Z) = rV^T^ ' ] = (z* + r«)3/2 ' (21) 

where the primes indicate the partial derivative with respect to I. Additionally, the unit 
normal vector reads 

n(r, z) = - 1 I Z ] . (22) 
V ; Vz' 2 + r' 2 \-r'j 

The differential equation ([TKjl is then solved by a forward relaxation scheme in time t 

r(t + At) \ = / r(t) \ _ Atfi fa z ) de fa z ^ (23) 
z(t + At) I \z{t) I 

where the steady-state solution d(r, z)/dt = is the solution of de(r, z) = we are looking 
for. In the numerical calculations we use a grid of 500 bins and an integration time step of 
At =0.001. The first and second derivatives are approximated using a symmetric two and 
three-step finite difference equation, respectively. Convergence is usually reached after 10 5 
time steps. The result is observed to be independent of the initial choice of r(z) at t = 0. 



Appendix B: Numerical solution of the PB equation 



Since we neglect mobile ions in our work, the PB equation reduces to Poisson's equation. 
It is solved on a two dimensional grid in cylindrical coordinates r and z with a finite difference 
method. The gradient and Laplacian are given then by V = (dr, dz) and A = d r + d r /r + d^, 
respectively. The first and second derivatives are approximated using symmetric two or 
three-step finite-difference equations. An explicit, forward time relaxation scheme is used to 
find the solution of Poisson's equation: 

(t + At; r) = *(t; f) - AtPB(^(t; f)). (24) 

In most cases we use a lattice spacing of Ar = Az = 0.4A on a n r x n z = 100 x 200 grid, 
and an integration time step At = 0.05. Convergence takes approximately 10 5 time steps. 
For the charged particles which are buried in the nanosolutes we use homogeneously charged 
spheres with a radius of 2A. Instead of a sharp transition for the dielectric boundary (|12p. 
we use a smoothing function for reasons of numerical stability: 

e(rO= r^YXl + (25) 
exp(Kd{r)) + 1 

where the absolute value of the length d{f) is given by the nearest distance to the boundary 
of the volume exclusion function v(r). d is defined to be positive when r e V and negative 
elsewhere. The inverse length k defines the width of the boundary region and in the limit 
k — > oo we recover the sharp transition (|12j) . We choose a value k > 3A 1 for which the 
solution of Poisson's equation becomes basically independent of the choice of k. An example 
for the dielectric boundary is shown in Fig. |H] for two partly dewetted nanosolutes of radius 
Ro = 15A at a distance sq = 7 A carrying a charge Q = 5e (system V in sec. III.C). 

In order to obtain the optimal shape function v m i n (f) the shape equation (fTfi|) has to be 
solved simultaneously with Poisson's equation when the solutes are charged. In practice, 
we first solve (|16J1 without any electrostatic contributions. In the second step, we solve 
Poisson's equation with the dielectric boundary (J25j) given by the volume exclusion function 
of the former solution. The result for the electric energy density is then plugged back in 
the shape equation in the third step. The last two steps are repeated until the solution for 
^min(^*) is fully converged. Since the results for r(z) excluding and including electrostatics 
are quite similar for our systems, full convergence takes usually only 6 to 7 repetitions of 



the described iteration steps. 
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e/(kJ/mol) 


cr/A AG sim /(kJ/mol) 


<Jbf/A -Rmin/A 


SPC 


0.65 


3.17 
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SPC/E 


0.65 


3.17 
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Na° 


0.2005 


2.85 


9.2(1) 


0.79 


2.32 


K° 


0.0061 


4.52 


23.7(5) 


0.76 


2.83 


Ca° 


0.6380 


3.17 


10.2(3) 


0.85 


2.80 


F° 


0.5538 


3.05 


9.7(2) 


0.85 


2.68 


Cl° 


0.5380 


3.75 


21(3) 


0.80 


3.30 


Br° 


0.4945 


3.83 


24(3) 


0.77 


3.35 
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Ne 


0.3156 


3.10 


11.41(0.05) 


0.84 


2.61 


Ar 


0.8176 


3.29 


8.68 (0.08) 


0.90 


2.96 


Kr 


0.9518 


3.42 


8.12 (0.1) 


0.91 


3.10 


Xe 


1.0710 


3.57 


7.65 (0.15) 


0.92 


3.27 














Me 


0.8941 


3.44 


10.96 (0.46) 


0.85 


3.11 


ethane 


— 




10.75 (0.50) 


0.87 




CH 3 


0.7503 


3.46 








propane 






13.81 (0.54) 


0.94 




butane 






14.69 (0.54) 


0.96 




CH 2 


0.5665 


3.52 








CH 3 


0.6900 


3.52 









TABLE I: Solute-water LJ parameters and solvation free energy AG s i m for neutral Lennard-Jones 
spheres from the SPC water simulations performed by Hummer et al£* and Pascheki^l <5 bf is the 
Tolman length best fit to AG s ; m (rounded to two digits after the decimal point). i? m i n is the 
resulting optimal radius excluded of solvent. Also shown are the values for the simple alkanes 
methane (Me), ethane, propane, and butane from the study of Ashbaugh et a/. 13 Simulation errors 
are given in parentheses. 



Atom 


AG smi / (kJ mol 


- 1 ) 5 b{ /k i? min /A 


Ne 


11.65 (0.05) 


0.88 2.60 


Ar 


8.83 (0.08) 


0.96 2.94 


Kr 


8.20 (0.1) 


0.98 3.09 


Xe 


7.58 (0.15) 


1.00 3.25 



TABLE II: Solvation free energies for neutral Lennard- Jones spheres in SPC/E water from the 
simulations of Pascheki^i 5m and R m i n are defined as in Tab. I. 



Ion 


q 




AG 


AG € 


Rmin/ A 


kJ mol -1 


kj mol - 1 


kJ mor 1 


Na + 


1 


-398 


-334 


-394 


1.83 


K+ 


1 


-271 


-246 


-282 


2.35 


Ca 2 + 


2 


-1306 


-1181 


-1364 


2.09 


F 


-1 


-580 


-274 


-630 


2.25 


cr 


-1 


-371 


-198 


-342 


2.97 


Br- 


-1 


-358 


-192 


-328 


3.02 



TABLE III: Solvation free energies for charged LJ spheres in SPC water from the simulations of 
Hummer et air 9 compared to the theoretical result AG. For 5 we use the best fits <5bf to the 
solvation of neutral spheres as shown in Tab. I. AG^ is the result when the dielectric boundary 
shift £ is applied, see text. 



System 


s/k 


vdW attraction 


Z 


W(s )/k B T dewetted 


I 


0.00 


no 





-57.6 


yes 


II 


0.75 


no 





-34.1 


yes 


III 


0.75 


yes 





-6.3 


yes 


IV 


0.75 


yes 


4 


-9.2 


yes 


V 


0.75 


yes 


5 


-5.1 


no 


VI 


0.75 


yes 


1 (oc) 


-1.3 


no 



TABLE IV: Studied systems for two alkane-assembled spherical solutes. W(sq) is the inter-solute 
pmf. If r(z = 0) 7^ the system is 'dewetted'. In system VI the solutes' charge is located off-center 
(oc) at the solute surface. 



FIG. 1: The particular solvation energy contributions AGi(R) with i = p,int,ne in Eq. (|18jl for 
one LJ sphere with Na° parameters given in Tab. I. The pressure term AG pr (thin solid line) 
with P = latm is basically zero on this scale. The interfacial term AG- m t(R) (dotted line) with 
7iv = 65m J/m 2 increases with radius R. The LJ term AG nc (R) is given by the dashed line. The sum 
of the three contribution gives the total AG(R) (solid line) with a minimum at R min = 2. 32 A for 
the uncharged sodium Na°. The inset shows the electrostatic contribution AG CS (R) (dot-dashed 
line) and the total AG(R) for the charged Na + with a minimum at R min = 1.83A. The best-fit 
Tolman length is <5bf = 0.79A. 



FIG. 2: Mean H(z) and Gaussian K{z) curvature and shape function r(z) (solid lines) for ethane. 
The canonical SAS (dashed line) from rolling a probe sphere with radius r p = 1.4A over the vdW 
surface (shaded region) is also shown. The vdW surface is defined by the solute-solute LJ-radius 
a ss /2 = 1.73AiSa 



FIG. 3: Mean H(z) and Gaussian K(z) curvatures and shape function r(z) for two alkane-assembled 
solutes of radius i?o = 15A (shaded region) at a distance sq = 8A for systems I- VI. The position 
of the charges Z = ±1 in VI are indicated by arrows. Curvatures are not shown for the 'wet' 
systems V and VI. 



FIG. 4: Top frame: theoretical pmfs for the systems I-III, and VI versus the solute distance sq. 
Bottom frame: corresponding mean forces. 



FIG. 5: Mean force /3FA between to nanosized solutes versus surface-to-surface distance so- 
The symbols denote the MD simulation results from Dzubiella et aZ. 29 ' 30 for neutral spheres with 
radius Rq = 12A (circles) and Rq = lOA(squares), and oppositely charged spheres with radius 
Rq = 10 A and charge Q = 2e (diamonds) and Q = 5e (triangles). The corresponding theoretical 
results using 5 = 0.9A are shown by solid lines; the range of the strong hydrophobic attraction 
decreases with decreasing radius and increasing charge. Dotted lines through the symbols are 
guides for the eye. The dashed line is the theory for Rq = 12A and 5 = 0.75A. 



FIG. 6: Distribution of the dielectric constant in space for two nanosolutes with Rq = 15A at 
a distance sq = 7 A carrying a charge Q = 5e (system V). The region between the spheres is 
dewetted. The distribution is scaled by ei = 78. 
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